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Abstract: Nonparametric or distribution-free charts can be useful in statisti- 
cal process control problems when there is limited or lack of knowledge about 
the underlying process distribution. In this paper, a phase II Shewhart-type 
chart is considered for location, based on reference data from phase I analysis 
and the well-known Mann- Whitney statistic. Control limits are computed us- 
ing Lugannani-Ricc-saddlepoint, Edgeworth, and other approximations along 
with Monte Carlo estimation. The derivations take account of estimation and 
the dependence from the use of a reference sample. An illustrative numeri- 
cal example is presented. The in-control performance of the proposed chart is 
shown to be much superior to the classical Shewhart X chart. Further com- 
parisons on the basis of some percentiles of the out-of-control conditional run 
length distribution and the unconditional out-of-control ARL show that the 
proposed chart is almost as good as the Shewhart X chart for the normal 
distribution, but is more powerful for a heavy-tailed distribution such as the 
Laplace, or for a skewed distribution such as the Gamma. Interactive soft- 
ware, enabling a complete implementation of the chart, is made available on a 
website. 



1. Introduction 

Control charts are most widely used in statistical process control (SPC) to detect 
changes in a production process. In conventional SPC, the pattern of chance causes 
is often assumed to follow the normal distribution. It is well recognized however that 
in many applications the underlying process distribution is not known sufficiently 
to assume normality (or any other parametric distribution) , so that statistical prop- 
erties of commonly used charts, designed to perform best under the assumed dis- 
tribution (such as normality), could be potentially (highly) affected. In situations 
like this, development and application of control charts that do not depend on nor- 
mality, or more generally, on any specific parametric distributional assumptions, 
seem highly desirable. Distribution-free or nonparametric control charts can serve 
this purpose. Chakraborti, Van der Laan and Bakir [4] (hereafter CVB) presented 
an extensive overview of the literature on univariate nonparametric control charts. 
In-control (stable) properties of these charts are completely determined (known) 
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and remain the same for all continuous distributions and hence their out-of-control 
behavior are more meaningful and comparable. 

In this paper we use the well-known Mann- Whitney test statistic as a charting 
statistic for detecting location shifts. The problem of monitoring the center or the 
location of a process is important in many applications. The location parameter 
could be the mean or the median or some percentile of the distribution; the latter 
two are often more attractive when the underlying process distribution is expected 
to be skewed. Among the available control charts for the mean of a process, the 
classical Shewhart X chart is the most popular because of its inherent simplicity 
and practical appeal. In some applications, the process distribution and/or the 
parameters are specified or can be assumed known. This is typically referred to as 
the standards known case (Montgomery [11], page 228). If on the other hand the 
parameters are unknown and are estimated from data, there is growing evidence in 
the recent literature that most standard charts, including the Shewhart X chart, 
behave quite poorly in terms of the false alarm rate and the average run length. 
We do not assume that the process parameters are specified or that the process 
distribution is known, instead, we assume that a reference sample is available from 
the in-control process from a phase I analysis. Once the control limits arc determined 
from the reference sample, monitoring of test samples is begun. This is referred to 
as a phase II application. 

There are some phase II nonparametric charts available in the literature; the 
reader is referred to CVB for many references and detailed accounts. The non- 
parametric charts considered by Chakraborti, Van der Laan and Van de Wiel [5] 
(hereafter CVV) are based on the precedence test. It is seen that the precedence 
charts arc good alternatives to the X chart in some situations. However, while the 
precedence charts are a step in the right direction, it is known that the nonpara- 
metric test underlying this chart, the precedence test, is neither the most powerful 
test (for location) nor the most commonly used nonparametric test in practice. 
With this motivation, we consider a chart based on the popular and more powerful 
Mann- Whitney [9] (hereafter MW) test which is equivalent to the perhaps more 
familiar Wilcoxon [15] rank-sum test. This is called the MW control chart. 

One might suspect that the distribution-free-ness of the MW chart might come 
at a "loss of power" with respect to parametric charts. However, remarkably, even 
when the underlying distributions are normal, the MW test is about 96% as ef- 
ficient (Gibbons and Chakraborti [7], pages 278-279) as (the most efficient) t-test 
for moderately large sample sizes, and yet, unlike the t-test, it does not require nor- 
mality to be valid. Park and Reynolds [12] realized the potential of nonparametric 
control charting and introduced a chart based on this statistic. They considered 
various properties of this chart when the reference sample size approaches infinity, 
which essentially amounts to assuming the standards known case. While this is 
important for theoretical purposes and gives some insight, such a chart does not 
appear to be very useful in practice where parameters and/or the underlying dis- 
tribution are unknown and need to be estimated from a moderate size phase I data 
set. In fact, it is crucial to develop and implement the MW chart in practice for 
small to moderate reference sample sizes since, as we show later in the paper, the 
MW chart can be especially useful in such cases. 

While the principles of this chart are simple, practical implementation of the 
chart, i.e., developing an efficient algorithm for computing the control limits, re- 
quires some effort. We provide software for practical use of the chart. Effectiveness 
of the chart is examined on the basis of several in-control and out-of-control per- 
formance criteria. We conclude with a discussion, including some topics for further 
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research. 

2. The MW control chart 

Suppose that a reference sample of size m, denoted by X = (X±, . . . , X m ), is avail- 
able from an in-control process and that Y = (Yi, . . . ,Y n ) denotes an arbitrary 
test sample of size n. The superscript h is used to denote the h th test sample, 
Y h = (Yi, . . . , Y£), h — 1,2, ... , when necessary for notational clarity; otherwise, 
the superscript is suppressed. Assume that the test samples are independent of each 
other and are all independent of the reference sample. The MW test is based on the 
total number of (X, Y) pairs where the Y observation is larger than the X. This is 
the statistic 

m n n 

(2.1) M XY = Y / Y1 r ( X * < Y i) = E >*i)+-" + I(Y 3 > X m )}, 

i=l j=X j=l 

where I(Xi < Yj) is the indicator function for the event {Xi < Yj}. Note that 
Mxy lies (attains values) between and mn and large values of Mxy indicate a 
positive shift, whereas small values indicate a negative shift. 

The proposed two-sided MW chart uses M XY as the charting statistic, which is 
Mxy for the h th test sample. The chart signals if 

Mxy < L mn or M XY > U mn , 

where L mn and U mn are the lower control limit (LCL) and the upper control limit 
(UCL), respectively The distribution of Mxy is known to be symmetric about 
mn/2 when the process is in-control, so it is reasonable to take L rnn = mn — U„ L n- 
We focus on two-sided charts, one-sided charts can be developed similarly. 

2.1. Design and implementation 

Implementation of the chart requires the control limits. Typically, in practice, 
the control limits arc determined for some specified in-control average run length 
(ARLq) value, say 370 or 500. If the successive charting statistics M XY , M| y , . . . 
corresponding to test sample 1,2,... were independent, then, as in the standards 
known case, the ARLq would be equal to the reciprocal of the false alarm rate, 
po = 2Pq(Mxy > Umn), where the subscript denotes the in-control case. So if the 
charting statistics were independent, the upper control limit U mn would be equal 
to the two-sided critical value for a MW test with size equal to 1/ARL - However, 
such critical values are not expected to be found in available tables for the MW 
test, since in a typical control charting application ARLq = 370, which means that 
the UCL is the upper critical value for a MW test with size 1/2(1/370) = 0.00135. 

Even if such critical values could be found, the main problem with their use is 
that the successive charting statistics M XY , M XY . . . are dependent, since the test 
samples are all compared to the same control limits derived from the same reference 
sample, and this dependence affects all operational and performance characteristics 
of the control chart, such as the false alarm rate, the ARL, etc. (see, for example, 
Quesenberry [13] and Chakraborti [3] for the Shcwhart X chart). It might be argued 
that for "large" amounts of reference data such dependence can be ignored. There 
are two problems with this argument, however. One, we would need to know the size 
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of the reference data, a-priori, that would support ignoring dependence, and two, 
we would have to wait much longer to gather that amount of data while process 
monitoring has to wait, costing time and money. The solution to this is the ability 
to calculate the control limits for any given (fixed) m and n and ARLq in a given 
situation. To this end we first develop an efficient method to calculate the ARL. 



2.2. Calculation of the ARL 



Let F and G represent the cdf of X and Y respectively, and suppose that F and G 
are continuous, so that "ties" between the X's and the V's, as well as within the 
X's and the Y's can be ignored theoretically. It is convenient to derive the ARL by 
conditioning on the reference sample, i.e. using the so-called conditioning method. 
To this end, observe that the probability of signal for any test sample, given the 
reference sample X = x, is 

(2.2) PG (x) = P G (M xY < mn - U mn ) + P G (M xY > U mn ). 

Let N denote the run length random variable for the chart. Given the reference sam- 
ple X = x, and that two arbitrary test samples Y h and Y l I) are independent, 
which implies the independence of M^ Y and M l xY . Hence, 

ARL = E(N) = E F [E G (N\X = x)} = #f(— ^) 

Pg{x) 



OO oo 



(2.3) 



/ ••• / — ^ dF{x x ) ■ ■ ■ dF(x m ) 
J J Pg{x) 

v(G( Xl ), . . . , G(x m )) dF{xi) ■ ■ ■ dF(x m ), 



-oo — oo 
oo oo 



say. The second equality in (2.3) follows from a property of expectation. The third 
equality follows since given X = x, N is geometrically distributed with parameter 
Pg{x). The fourth equality is obtained by writing l/po(x) as a function of G and 
x\, . . .,x m , say, v(G(xi), . . . ,G(x m )) = u{P G {Y :j < Xl ), . . . ,Pc{Yj < x m )), where 
v is some function. This can be done since pq(x) is a sum of probabilities like 
Pg(M x y = u )i which, in turn is a sum of probabilities over all configurations 
of the x's and y's for which M x y equals u. Naturally, the probability of such a 
configuration only depends on G(x±), . . . , G(x m ). 

The (unconditional) ARL of the chart is the mean (expectation) of the distribu- 
tion of the random variable Eg(N\X), which is the conditional average run length, 
given the random reference sample. Percentiles of the distribution of Eg(N\X) (and 
not just the mean) are useful to study and characterize control chart performance 
when parameters arc estimated, and we develop efficient algorithms to compute 
these. First however, we focus on the mean of the conditional distribution, that is 
the unconditional ARL given in (2.3), when the process is in-control. 

In the in-control situation the X's and the y's come from the same distribution 
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F = G. Therefore, we may assume w.l.o.g. that F = U[0, 1]: 

oo oo 

ARL = I ■ ■ I v{F{ Xl ),...,F{x m )) dF( Xl )---dF(x m ) 



— oo — oo 

1 



(2.4) = / ■ ■ • / v(u\, . . . , u m ) dm ■ ■ ■ du % 

1 





1 1 







pu(u) 



dui ■ ■ ■ du m , 



where pu( u ) = Pu(MuY < Tnn — U mn ) + Pjj{M u y > U mn ) is the conditional 
probability of a signal at any test sample, given the reference sample u, when 
the process is in-control. The subscript U is used to denote that in the in-control 
case both the reference and the test samples can be thought of coming from the 
same distribution, the C/[0, 1] distribution, which shows that the in-control ARL 
of the proposed chart does not depend on F. The same argument can be used to 
show that the in-control run-length distribution does not depend on F and hence 
the proposed chart is distribution-free. We emphasize that for the in-control case, 
without any loss of generality, the common (F = G) distribution can be assumed to 
be the U[0, 1] distribution by virtue of the probability integral transform (see, for 
example, Gibbons and Chakraborti [7]). This simplifies calculations considerably. 

We need to calculate (2.4) to implement the chart and (2.3) to evaluate chart 
performance. For both of these objectives there are two problems. First, an explicit 
formula for pg{x) (or pu(u)), is not known, which prevents a direct computation. 
Second, we have an m-dimensional integration in both (2.3) and (2.4). Our approach 
is to calculate pc(x) (and pu(u)) (exactly or approximately) using a fast algorithm, 
and then use that to approximate the integral in both (2.3) and (2.4) with a Monte 
Carlo estimate to get estimates 

1 - 1 

(2.5) ARL G n-"£^^ 
and 



1 K 

(2.6) ARL 



f Pu{ui) ' 

where x t = (xn, . . . , x, m ) is the i th Monte Carlo sample, i = 1, . . . , K, of which 
each component is drawn from some specified F for the ARLq , and K denotes the 
number of Monte Carlo samples used. Similarly, for the in-control situation, u t is 
a Monte Carlo sample from U[0, 1]. For an accurate approximation, K needs to be 
sufficiently large and therefore a fast method of computing the signal probability 
Pg{x) (and pu{u)) for an arbitrary reference sample x is essential for the practical 
use of the approximation. 

The ARL calculations proceed in two steps. The first step is to find a fast and 
efficient method to compute the signal probability. We detail the procedure for fast 
computation of Pg{x)', calculation of Pu(u) is similar. 
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2.3. Fast computation of signal probability 

From (2.1) it is seen that the MW statistic is the sum Y^j=i where Cj repre- 
sents the number of X's that are less than an Yj. From (2.2) it follows that the 
calculation of pc(x) essentially requires calculation of the upper-tail probability 
Pq{M x y > U mn ), say, and this in turn requires (i) efficient enumeration of all n- 
tuples {C\, . . . , C„} for which the sum is greater than U mn and (ii) summation of 
the probabilities for such tuples. 

Note that P(Cj = I) is equal to P(Xm < Yj < Xn + W), where Xm denotes 
the I th ordered observation in the reference sample for I = 1, . . . , m with X/ ) = 
— oo and Xr m +i) = oo, say. Given the reference sample, the last probability is 
simply P(xm < Yj < xn+u), which is denoted by a;, I = 0, . . . , m. Also, given the 
reference sample, note that the random variables Cj are i.i.d. Hence the conditional 
probability generating function (pgf) of Cj is 

rn m 

(2.7) Ht(z) =^P(Cj =l)z l =5^ojz'. 

1=0 1=0 

Again, since given the reference sample the Cj are i.i.d., the conditional pgf of M x y 
(the sum of the Cj), is simply the product of the pgf's in (2.7) 

mn rn 

(2.8) H 2 (z) = £p(M x y = j)^ = (E^')" 

3=0 3=0 

In principle, Pg(M x y > U mn ) can be calculated by expanding the power in (2.8) 
and collecting the coefficients of all terms with degree greater than U mn . However, 
for moderate to large m (say, m > 100) and n not very small (say, n > 5) this takes 
a considerable amount of computing time, especially since the procedure has to be 
repeated K (a large number) times, once for each Monte Carlo sample. Alternative, 
faster, methods such as "a branch-and-bound" algorithm, based on Mehta et al. [10] 
can be used to shorten the intermediate expressions that result from expanding (2.8) 
and saves considerable time. 

However, even with the branch-and-bound algorithm, m and n might be just too 
large to allow for exact computations and hence a good approximation to the control 
limits may be necessary. For example, we may apply the central limit theorem for the 
sum of i.i.d. random variables to M xY = Y^j—i Cj to get a normal approximation 
to Pg(M x y > U mn ) but in our context, U mn is typically far in the upper tail 
of the distribution of M x y and n is usually not very large, and so the normal 
approximation is not very accurate. Instead, we find the Lugannani-Rice formula 
(hereafter LR- formula; see Jensen [8] , page 74) for the upper-tail probability for the 
mean of i.i.d. discrete random variables (which is a "saddlepoint" approximation 
formula) to be particularly useful. This formula is known to be more accurate than 
the normal approximation in the tails of a distribution and is based on the cumulant 
generating function of Cj, which is obtained from the pgf in (2.7): k(t) = log[i?i(e*)]. 
Let m(t) and cr 2 (t) denote the first and the second derivative of k(t), respectively. 
Furthermore, let u = (U mn + l)/n and M x y = M x y/n. The saddlepoint 7 is the 
solution to the equation m(t) = u. Using (3.3.17) in Jensen (1995, page 79) we 
obtain 

P G (M xY >U mn ) = P G (M xY > U mn /n) = P G (M xY > u) 

w l-$(r) + 0(r) (---), 
A r 
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where 

A = n 1/2 (l - e 7 )cr(7), r = (sgnj){2n(ju - fc( 7 ))} 1/2 , 

$(.) and </>(.) arc, respectively, the cdf and the pdf of the standard normal distribu- 
tion. Using (2.9), we can efficiently approximate the signal probability pg(%) given 
in (2.2). 

2. 4- Monte Carlo estimation of ARL and error control 

From formulas (2.5) and (2.6) we observe that the computation of Pg(%) is repeated 
many times to obtain a Monte Carlo approximation of the ARL. The question here 
is regarding K, the number of samples so that the Monte Carlo error is acceptably 
small. Since, for the purpose of ARL computation, the reference samples are drawn 
independently from G (or U[Q, 1] in case of ARLq), the Monte Carlo standard error 
is estimated by s mc = s(ARLg(X)) /y/~K ', where s(.) denotes the sample standard 
deviation computed from K simulated reference samples. Then, we may choose the 
smallest K such that 

(2.10) Smc = s{ARL G {X)) / VK < D, 

where D is either a specified number or a percentage of the current estimate ARLq. 
We can start with say K = 100, increase K, compute s mc , and repeat the process 
until the specification s mc < D is met. Use of formula (2.10) provides a way to 
obtain an efficient and a reasonably accurate approximation to ARLq- Asymptotic 
probabilistic control of the accuracy may be achieved by using the normal distribu- 
tion for ARLq. For example, one could set D such that the probability that ARLq 
deviates more than C units from the real mean is smaller than 0.05. Here, C could 
be a small percentage of the current estimate. Next, we discuss the approximation 
of ARL in more detail. 

Approximation of ARLq 

Three methods of approximating ARLq have been introduced so far, each based on 
a different method to compute or approximate pu(u). To summarize, they are: 

1. Exact (EX): Monte Carlo simulation using (2.6), with pu{x) computed exactly 
using formula (2.8) 

2. LR-formula (LR): Monte Carlo simulation using (2.6), with pu(u) computed 
approximately using formula (2.9) 

3. Normal (NO): Monte Carlo simulation using (2.6), with pu(u) computed from 
a normal approximation 

We compare these methods on the basis of accuracy and speed. While setting 
a value of D, we observed that K was usually under 1000, and the maximum 
Monte Carlo error was 2% of the target ARLq = 500, that is, equal to 10. The 
value of K was set to 1000 and kept unchanged in these computations to get a 
fair comparison of the computing times. Moreover, we introduce two alternative 
methods to calculate ARLq: 

4. Fixed reference sample (FR): Fix reference sample to q = (l/(m + 1), . . . , 
m/(m+ 1)) and approximate ARLq by l/pu(q) 

5. l/(falsc alarm rate) (FA): approximate ARLq by the reciprocal of the false 
alarm rate: 1 /(2Pq(Mxy > U mn )). 
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Table 1 

ARLq approximations and computing times 



EX LR NO FR FA 







ARL 


time 


ARL 


time 


ARL 


time 


ARL 


time 


ARL 


Time 


m 


n 




(sec.) 




(sec.) 




(sec.) 




(sec.) 




(sec.) 


50 


5 


486 


54 


506 


36 


307 


1.0 


403 


0.05 


247 


0.01 




10 


504 


395 


505 


34 


327 


1.0 


524 


0.05 


226 


0.01 




25 


488 


4850 


491 


31 


425 


1.2 


694 


0.05 


119 


0.01 


100 


5 


496 


220 


505 


48 


219 


1.2 


478 


0.05 


353 


0.01 




10 


505 


1920 


506 


47 


339 


1.3 


531 


0.05 


332 


0.01 




25 


** 


26168 


503 


18 


422 


1.3 


683 


0.06 


233 


0.01 


500 


5 


491 


10633 


496 


207 


226 


1.2 


492 


0.20 


445 


0.01 




10 


** 


73516 


513 


179 


367 


1.7 


537 


0.21 


484 


0.01 




25 


** 


7.59*10 5 


494 


176 


445 


1.6 


578 


0.29 


450 


0.01 


1000 


5 


** 


31766 


500 


356 


235 


2.1 


513 


0.48 


171 


0.01 




10 


** 


3.42*10 5 


499 


373 


355 


2.4 


516 


0.49 


488 


0.01 




25 


** 


3.15*10 6 


500 


348 


442 


1.7 


548 


0.63 


482 


0.01 


2000 


5 


** 


1.71*10 5 


503 


713 


234 


2.1 


506 


0.67 


474 


0.01 




10 


** 


1.44*106 


504 


659 


354 


1.9 


513 


0.71 


499 


0.01 




25 


** 


1.29*10 7 


509 


676 


446 


2.1 


531 


1.41 


497 


0.01 



**ARLo could not reliably be estimated within reasonable time; computing time for K = 1000 is 
obtained by multiplying computing time for K = 1 by 1000 (sampling algorithm is linear in K). 



Let us explain approximations 4 and 5. When m is large, the empirical cdf F m {x) 
converges to F{x) (which is the cdf of the U[0, 1] distribution) and hence for large to 
we may approximate the i th reference sample observation by the i / (m+l) th quantile 
of the U[0, 1] distribution, qt = i/(m + 1), i = 1, . . . , to. Thus, we approximate the 
ARLq, for large m, by l/pu(q), where q = (q\, . . . , q m ). This is method 4. The major 
benefit with this approximation is that we need to compute pu (u) only once (namely 
at u = q) instead of K times as needed in methods 1 through 3 for each of K Monte 
Carlo reference samples. Finally, another quick approximation for the ARLq is given 
by the inverse of the false alarm rate, 1 /(2P (Mxy > U mn )). In a setting where 
runs are truncated at a finite point T, the latter approximation was proven to be 
unbiased when m approaches infinity in Park and Reynolds [12] . This approximation 
is method 5. Chakraborti [3] showed that for the Shewhart X chart, 1/FAR is a 
lower bound to ARLq and noted that this bound can serve as a "quick and dirty" 
approximation to the ARLq for moderate to large values of to. When applying 
method 5, we used formula (11) in Fix and Hodges [6] to compute Pq(Mxy > U mn ), 
based on an Edgeworth approximation, which significantly improves the normal 
approximation by including moments of order higher than 2. 

Tabic 1 displays the estimated ARLq values computed by the five methods for 
fifteen combinations of to and n. Chart constants were determined (using the al- 
gorithm discussed in the next section) such that ARLq f=a 500, when applying the 
exact formula (2.8) or the best approximation, the LR-formula, when exact compu- 
tations were too time-consuming. Therefore, the closer an ARLq value is to 500, the 
better the approximation. The table also shows the computing times on a 1.7GHz 
Pentium PC with 128MB of internal RAM. 

Several observations can be made from Table 1. First we see that the "gold 
standards" or the exact computations are very time-consuming for most values 
of to and n. However, when they can be found, they would naturally form the 
basis of our comparisons of various approximations. Second wc sec that the normal 
approximation is not very accurate, but the LR-approximation is, particularly for 
to < 100, n < 10 and for to = 500, n = 5. Since the LR-approximation is known 
to become more accurate when the sample sizes increase, one may safely apply 
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the LR-formula also when to > 50 and n > 5 in order to implement the proposed 
chart. It may be noted that when to increases, the computing times with the LR- 
formula also increase, although, by far, not as dramatically as the times for the 
exact computations. This suggests that in practice (for finding the chart constants, 
to be discussed next) there is still room for an alternative, quick approximation of 
ARLq, particularly for large values of to. Compared to the LR-formula, we observe 
that both the "fast" approximations (methods 4 and 5) are quite good for to > 1000 
and that the fixed reference sample approximation (method 4) performs somewhat 
better for relatively small values of n (n = 5,10) than for n = 25. 

To summarize, the best method of calculating the ARLq is the exact EX method 
if it is computationally feasible, otherwise, the best approximation is the LR meth- 
od. In practice, we recommend using the LR-method, because it is both fast and 
accurate. If the reference sample is very large, say to > 1000, one of the two faster 
approximations, either FR or FA, can be used. 

2. 5. Determination of chart constants 

Since we can now calculate the ARLq corresponding to a given value of U mn effi- 
ciently and accurately, we can use an iterative procedure based on linear interpola- 
tion to find the control limit for a pre-specified ARLq value, say 500. In principle, we 
use the LR-approximation for the computation of ARLq. However, we have observed 
that this approximation is still somewhat time-consuming, so we want to minimize 
the number of iterations for which the LR approximation is used. To this end a 
good starting value of U mn is needed and this is where the fast approximations FR 
and FA are very useful. Starting with the FA approximation, we simply equate the 
inverse of the false alarm rate to 500, which means solving 1/(2* FH(u)) — 500 for 
u in order to get an initial guess for U mn . Since the FR approximation is somewhat 
better than FA for to < 500, we use this initial guess with the FR approximation to 
refine the Fix-Hodges approximation when to < 500. The resulting approximation 
for the UCL is then used as an initial guess for the linear interpolation method 
with the LR-formula. We do not detail the search procedure here, but illustrate it 
with an example. Suppose that to = 375 and n = 7, and we want to find the chart 
constants such that ARL ~ 400. Suppose we allow a deviation of 2% (which is 0.02 
* 400 = 8) maximally, hence the search procedure would stop and yield the desired 
control limit U mn , when 392 < ARLq < 408. Moreover, suppose we stipulate that 
the Monte Carlo standard error s mc be smaller than 1.5% of the current estimate 
of ARLq; from inequality (2.10) wc observe that this requirement determines the 
number of Monte Carlo samples K per iteration, when setting D— 0.015 * 400 = 6. 
The output from our program (written in Mathcmatica; sec Software section later) 
is shown in Table 2. As can be seen in Table 2, six iterations (numbered 1 through 6) 
have been executed in approximately 140 seconds. The first three of these hardly 
take any computing time, because the fast FA and FR approximations were used. 
For each iteration, the values of the UCL and the LCL and the corresponding ARLq 
are shown. Under the LR method, the program also calculates the 5 th percentile 
of the conditional in-control ARL distribution and the Monte Carlo standard error 
(smc). Note that the first step with the LR method, step 4, uses the chart constants 
of step 3, which is our best guess from the fast approximations, as initial values. 
Also, note that for the LR method, the first two iterations (4 and 5) produce ARLo 
values below and above the target value 400, so that linear interpolation begins 
at the third iteration and the new UCL is found using the two previous UCUs 
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Finding control limits for m = 375, n =7 and target ARLq = 400 

FA: l/(false alarm rate) approximation 

1. ucl=2146 lcl=479 ARL0=400 

FR: Fixed reference sample approximation 

2. ucl=2146 lcl=479 ARL0=446.761 

3. ucl=2136 lcl=489 ARL0=386.729 
LR approximation 

4. ucl=2136 lcl=489 ARL0=380.059 smc=5. 69018 5% perc=238.407 K=402 

5. ucl=2146 lcl=479 ARL0=438.11 smc=6.5647 5% perc=287.53 K=319 

6. ucl=2139 lcl=486 ARL0=394.496 smc=5.91419 5% perc=252.778 K=315 
139.962 Second 



Table 3 

Lower and upper MW control chart limits for selected values of m, n and ARLo 







ARL 


= 370 


ARL = 


500 


m 


n 


Lmn 


Umn 


Lmn 


Umn 


50 


5 


35 


215 


33 


217 




10 


115 


385 


111 


389 




25 


400 


850 


393 


857 


100 


5 


69 


431 


65 


435 




10 


231 


769 


224 


776 




25 


805 


1695 


793 


1707 


500 


5 


348 


2152 


328 


2172 




10 


1170 


3830 


1128 


3872 




25 


4081 


8419 


4016 


8484 


1000 


5 


698 


4302 


653 


4347 




10 


2344 


7656 


2268 


7732 




25 


8169 


16831 


8058 


16942 


2000 


5 


1397 


8603 


1309 


8691 




10 


4682 


15318 


4540 


15460 




25 


16392 


33608 


16145 


33855 



and their corresponding ARLo values. Thus, the final chart constants are found at 
iteration 6, U mn — 2139 and hence L mn — 486, with attained ARLo— 394.5. The 
5 th percentile of the conditional in-control ARL distribution, at this iteration, is 
found to be 252.78. So using the MW chart with UCL = 2139 and LCL = 486, 
the unconditional ARLo = 394.5 implies that when the process is in-control, on an 
average, a false alarm is expected every 395 samples. The 5 th conditional percentile 
implies that for 95% of all reference samples (that could have possibly been taken 
from the in-control process), the average run length is at least 253. 

Table 3 shows chart constants computed with this algorithm for a number of 
combinations of m and n. All cases are two-sided and ARLo equals either 370 or 
500. 

2.6. Numerical example 

Table 5.1 in Montgomery [11] gives a set of data on the inside diameters of piston 
rings manufactured by a forging process. Twenty-five samples, each of size five, were 
collected when the process was thought to be in-control. The traditional Shcwhart 
X and R charts provide no indication of an out-of-control condition, so these "trial" 
limits were adopted for use in on-line process control. 

For the proposed MW chart with m = 125, n — 5 and ARLo = 400, we find 
the upper control limit U mn = 540 and hence the lower control limit L mn = 85. 
Having found the control limits, prospective process monitoring in phase II begins. 
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2 3 4 5 6 7 8 9 10 11 12 13 14 15 

Fig 1. MW Chart for the Piston-ring data. 

Montgomery also gives (Table 5.2) fifteen additional samples from the piston-ring 
manufacturing process. 

These "test samples" lead to fifteen MW statistics calculated using Minitab: 
429.0, 333.0, 142.5, 370.5, 241.5, 410.5, 393.0, 240.5, 471.0, 486.0, 340.5, 561.0, 
575.5, 601.5 and 484.5. Comparing each statistic with the control limits, all but 
three of the test groups, 12, 13 and 14 are declared to be in-control. The control 
chart is shown in Figure 1. 

The conclusion from this chart is that the medians of test groups 12, 13 and 14 
have shifted to the right in comparison with the median of the in-control distribu- 
tion, assuming that G is a location shift of F. It may be noted that the Shewhart 
X chart shown in Montgomery [11] led to the same conclusion with respect to the 
means. Of course, the advantage with the MW chart is that it is distribution-free, 
so that regardless of the underlying distribution, the in-control ARL of the chart is 
roughly equal to 400 and there is no need to worry about (non-) normality, as one 
must for the X chart. To see how the MW chart compares with other available non- 
parametric charts, we calculated the distribution-free precedence chart (see CVV) 
for this data. We found LCL = 73.982 and UCL = 74.017 for the precedence chart, 
with an attained ARL sa 414.0. Consequently, the precedence chart declares the 
12*^ and the 14" 1 groups to be out of control but not the 13 th group, unlike both 
the MW and the Shewhart chart. 

Comparison with the Shewhart chart 

The performance of a control chart is usually assessed in terms of its run length dis- 
tribution and certain associated characteristics, such as the ARL. While following 
the general norm in the literature we examine the ARL, given the skewed nature of 
the run length distribution, we also consider two other criteria for evaluating and 
comparing the performance of the MW chart and its parametric competitor, the 
Shewhart X chart, in terms of some percentiles of the conditional run length dis- 
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tribution. We believe these criteria provide additional useful information regarding 
chart performance, with estimated parameters. 

To ensure a fair comparison between the MW and the Shewhart A chart, first, 
the X chart is used for the case when both the mean and the variance are unknown, 
with parameters estimated from the reference sample. Second, the charts are both 
designed to have the same specified ARLq. Note that for the X chart the non- 
robustness of the ARLq with respect to non-normal in-control distributions is a 
major concern. This has been recognized as a problem elsewhere (see e.g. CVV) 
and is perhaps one of the important reasons for considering a nonparametric chart 
in practice. 

Although both the conditional and the unconditional distributions provide im- 
portant information regarding the performance of a control chart, we argue that 
the conditional distribution might be preferred from a practical point of view since 
the unconditional distribution is what results after "averaging" over all possible 
reference samples. Users would most likely not have the benefit of averaging in a 
particular application. Also, since the distribution is skewed, the percentiles and not 
the average are better measures of performance and thus the standard deviation 
appears to be a less suitable measure of variability. First we consider the in-control 
case. 

2. 7. In-control performance 

For the in-control situation, a lower order percentile (say the 5 th ) is more useful 
and relatively large values of this percentile are desirable (in the same spirit that 
the ARLq of a chart be large), since that would lead to a smaller probability of a 
false alarm. We also show the estimated standard deviations to give an indication 
of the variability of ARLq(X), since this appears to be the current norm. For 
completeness, the 95 th percentile is also given. The two percentiles can also provide 
an indication of the variability in the conditional distribution. 

The proposed MW chart for location for is compared to the Shewhart X chart 
with estimated parameters. There may be some interest in comparing against other 
control charts and we comment on this aspect later. To ensure a fair comparison, 
chart constants were determined such that the ARLq approximately equals 500 for 
both charts. We kept the test sample size constant, n = 5, and used several values 
for the reference sample size m. Both samples were drawn from a Normal(0,l) 
distribution. The number of simulations, K, was set to 1000. The results are shown 
in Table 4. 

To illustrate, suppose m = 750. While applying the MW chart with U mn = 3258, 
from Table 4 we know that 95% of the in-control ARUs (for a large number of ref- 
erence samples taken from the in-control process) will be at least 360. This provides 
useful performance information in addition to saying that the unconditional ARLq 
is 500. For the X chart with m = 750 and a control chart constant of 3.089 (this 
guarantees ARLq = 500 when parameters are both unknown), the 5 th percentile is 
314 so that 95% of the in-control ARUs are at least 314. Since the in-control 5 th 
percentiles for the MW chart are considerably larger than those of the Shewhart 
chart for all m, we conclude that the in-control performance of the MW chart is su- 
perior to that of the Shewhart chart with estimated limits, particularly for m < 150. 
Thus the MW chart is more useful in applications where a large amount of refer- 
ence data might not be available. The uniformly smaller standard deviations for the 
MW chart doubly confirm its superiority. Note also that as m increases, both the 



168 



S. Chakraborti and M. A. van de Wiel 



Table 4 

Fifth and ninety-fifth percentiles and standard deviations of the conditional in- control 
distribution of ARLq(X); All cases: n = 5 and ARLq = 500 



m 


Upper 
control 
limit 
MW 


5 th 
perc. 
MW 


95 th 
perc. 
MW 


St. 
Dev. 
MW 


Upper 
control 
limit 
Shewhart 


5 th 
perc. 
Shewhart 


95 th 
perc. 
Shewhart 


St. 
Dev. 
Shewhart 


50 


217 


97 


1292 


553 


3.01996 


49 


1619 


854 


75 


326 


146 


1219 


461 


3.05156 


87 


1379 


645 


100 


435 


182 


1146 


358 


3.06535 


112 


1290 


463 


150 


654 


251 


1090 


315 


3.07715 


154 


1197 


377 


300 


1304 


284 


845 


197 


3.08607 


232 


927 


235 


500 


2172 


322 


700 


140 


3.08848 


270 


828 


174 


750 


3258 


360 


677 


107 


3.08935 


314 


765 


140 


1000 


4347 


379 


674 


83 


3.08969 


338 


721 


121 


1500 


6520 


409 


642 


71 


3.08996 


367 


678 


97 


2000 


8691 


420 


629 


55 


3.09007 


376 


651 


84 



5 th as the 95 th percentiles approach the mean of the corresponding unconditional 
distribution, the ARLq, which is set at 500. 

2.8. Out-of -control performance 

The distribution-free property (and the resulting robustness of the ARLq) is an im- 
portant asset of the proposed chart, but what about its out-of-control performance? 
We address this issue here. Since this is a chart for location, our interest is in the 
shift (location) alternative G(x) = F(x — 8), where 6 is the unknown shift parame- 
ter. To study the effects of using a reference sample, again, wc use conditioning and 
study the distribution of the conditional ARL. Let ARLs(X) = E F ( x _g- ) (N\X), 
denote ARL of the run-length distribution, given the reference sample X , when the 
process distribution F has shifted by an amount 8. We examine the out-of-control 
performance of both the MW chart and the Shewhart chart in terms of ARLs(X) 
next. 

We also provide a more traditional chart comparison by examining the out-of- 
control unconditional ARL for a specified distribution and shift, ARL^. Naturally, 
small values of ARL^ are desirable. As in the case of the in-control situation, we 
also examine a percentile of the out-of-control distribution. However, in the out- 
of-control case it makes sense to focus on a higher order percentile, say the 95"* 
percentile. Denoting this by (?o.95, relatively smaller values of (?o.95 are desirable for 
a preferred chart, since the probability of a signal is desired to be higher in the out- 
of-control case. For a given value of go. 95, users can be 95% confident that for their 
own specific reference sample, the out-of-control ARL is smaller than that value. 
The two performance measures, namely the ARL^ and the go. 95, are examined for 
three distributions: Normal, Laplace and Gamma(2,2). The Laplace distribution is 
normal-like but with heavier tails, which results in higher probabilities of extreme 
values. The Gamma(2,2) distribution is skewed and is often used in the SPC lit- 
erature. We apply two-sided charts to the Normal and the Laplace distributions 
and a one-sided chart with an upper control limit to the Gamma(2,2) distribution. 
The test sample size n is 5 and the reference sample size m is 100. Control limits 
for both the MW chart and the Shewhart X chart with estimated parameters are 
determined such that ARLq = 500. Using these limits ARL^(A) and the 95 th per- 
centile of the distribution of ARL^(A) are computed for several values of 8, which 
is given in units of the standard deviation. Figures 2 through 4 show the results. 
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Fig 4. Performance for MW chart and Shewhart chart under Gamma(2,2) shift alternatives. 
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For the set of points, triangles and diamonds, observe that for the normal distri- 
bution the 95 th percentiles for the Shewhart X chart are all smaller than those for 
the MW chart. Thus, as one might expect, the X chart is more effective in detect- 
ing shifts than the MW chart in case of the normal alternative. However, note that 
the differences between the percentiles are small at all shifts (the largest difference 
is around 15) and the difference appears to vanish for shifts greater than 1. The 
same pattern holds for the two ARL's. On the other hand, Figure 3, for the Laplace 
distribution, shows that the MW chart is clearly better than the Shewhart chart 
for all shifts, large and small. For the Gamma(2,2) distribution, in Figure 4, again, 
we see that the MW chart is better in detecting shifts, although the difference in 
performance is not as dramatic as in the case of the Laplace distribution. These 
calculations were repeated for m = 500 observations; the results were very similar 
and are therefore omitted here. We conclude that the nonparametric MW chart 
follows the well-known results for the MW test statistic: it is nearly as effective as 
the Shewhart X chart under normality, but is more effective under heavy tailed 
and skewed distributions. Also, note that performance of the MW chart in the case 
of the Laplace distribution makes it potentially useful when outliers in the data are 
not uncommon. 

3. Discussions and further topics 

Comparison with CUSUM and EWMA charts 

Shewhart charts are known to be very good for moderate to large shifts. These 
charts do not require tuning parameters like those needed by the CUSUM and the 
EWMA charts for a specified shift, but aim for global performance. Thus although 
one can design a CUSUM or an EWMA chart to perform better by focusing on 
cither small, medium or large shifts, a priori, the same path could lead to a worse 
performance for other shifts. Moreover, as has been noted in the literature Quc- 
senberry [14], the very nature and principle of the "averages-type" charts (such 
as CUSUM and EWMA) is different from that of Shewhart charts. In addition 
to the problem of having different in-control run-length distributions that renders 
the stable properties of these charts to be quite different (and hence out-of-control 
assessments less meaningful), the averages charts are more powerful in detecting 
specific types of shifts (sustained monotonic) . 

Nevertheless, we made an attempt to compare the MW chart with EWMA and 
CUSUM on the basis of in-control robustness in terms of misspecifications of the 
shape of the distribution and the variance. Rather than showing all the results, we 
summarize the findings here. For fast detection of large shifts (say 1 and larger), 
the MW chart is very good: it is powerful and maximally robust against misspeci- 
fications of shape and variance. The latter is far from being true for CUSUM and 
EWMA designed for detecting large shifts: true ARL could easily be twice as large 
or small as the target ARLq in case of skewness, increased variance or heavy-tails. 
For small shifts, the situation is more delicate: the EWMA is a strong competitor, 
since it is often (but not always; for example, not for medium to large shifts in the 
normal) more powerful than the MW chart and quite robust against misspecifica- 
tions of the shape. This type of robustness of the EWMA chart was shown earlier 
(see e.g. Borror et al. [2]). However, the EWMA is seen to be not robust against a 
misspecification of the variance. For example, a small increase of the variance of a 
normal (from 1 to 1.1) lowered the in-control ARL from 500 to 215 for an EWMA 
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designed for small shifts. The overall conclusion is that especially when the in- 
control variance can be estimated with limited accuracy, the nonparametric charts 
in general, and the MW chart in particular, is extremely useful in practice, because 
they do not require knowledge of the underlying variance. Similar conclusions have 
been drawn in Amin ct al. [I]. 

Alternative chart design criteria 

Table 4 suggests that with estimated parameters, especially for small values of 
77i, it may be useful to use a lower order (say the 5 th ) percentile of the conditional 
distribution as a chart design criterion rather than the mean, i.e., the unconditional 
mean ARLq, The design requirement would be that the percentile be at least equal 
to some specified large number, such as 300. In some cases one might want to avoid 
very short in-control runs in the future, which suggests using a lower percentile 
of the in-control distribution of N and not that of the conditional distribution of 
E(N\X). The probability Pq(N < n) can be computed using similar methods as 
for computation of ARLq, Then, for example, if one wants to avoid in-control runs 
smaller than 100 with a high probability, say 0.90, we can find the control chart 
limits by solving Pq(N < 100) = 1 — 0.90 = 0.10, again using a search algorithm. 
To facilitate use of both of these chart design criteria we have implemented these 
in the software that we provide with this paper. 

Individual's chart 

There is considerable interest in nonparametric charts for individual observations 
(77 = 1) since in this case the notion of approximate normality via the central 
limit theorem is not applicable. For n = 1, formula (2.7) allows for fast exact 
computations. The software can be used to set up such a chart. Because of the 
natural interest in nonparametric individual's charts, a detailed treatment of this 
topic will be given in a future paper. 

Monitoring dispersion 

The dispersion or spread need not be monitored while using a nonparametric control 
chart (such as the MW chart) under the location model. This has been cited as an 
advantage of the nonparametric sign charts by some authors. However, monitoring 
the spread is an interesting practical problem in a "location-scale" model and we 
see the possibility of designing a chart for scale (along with that for the location) 
based on some nonparametric test for scale, This topic will be considered in a future 
paper. 

Appendix: Software 

In order to support practical implementation of the methods presented in this 
paper two types of software related resources are provided. First, a Mathematica 4.2 
Wolfram [16] notebook is available to calculate (i) the in-control ARL computations 
any of the five methods (ii) the out-of-control performance computations, (iii) the 
control chart constants and (iv) to plot the MW-control chart for a user-specified 
data set. Second, we created a website that enables anyone to apply the proposed 
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methodology. The MW control chart limits can be found for the sample sizes at hand 
for a specified target ARLq (or a desired (gth) percentile of ARLq(X)). Moreover, 
the website allows users to import their own data set and have the MW chart 
drawn. The site can be reached via www.win.tue.nl/~markvdw. The Mathcmatica 
notebook is available from the same site; it contains more procedures and allows 
for more flexible input. User instructions are available in the notebook and at the 
website. 
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